program shallow_water_main

    use common_info
    use space_manager
    use time_manager
    use mesh_manager
    use dynamic_quantities
    use dist_recon
    use constant
    use type_define
    use diagnostics
    use sample_data
    use message
    use run_manager

    implicit none

    character(20) :: namelist_file

    real, allocatable :: height_gradient(:,:)
    real, allocatable :: divergence(:,:)

    ! height field
    character :: height_name = "h" ! ---> This may need to be redesigned
    integer :: h_index
    real, pointer :: h(:)
    real, allocatable :: h_tendency(:)

    ! wind field
    character :: wind_name(2) = (/"u","v"/) ! ---> This may need to be redesigned
    integer, allocatable :: u_index(:)
    type(real_list), allocatable :: u_tendency(:)

    integer :: i, j, k, l

    ! parse command arguments
    call get_command_argument(1, namelist_file)
    if(namelist_file == "") then
        call message_show(ERROR_MSG, &
            "shallow_water_main",    &
            "shallow_water_main",    &
            "A valid namelist file is needed")
        call run_manager_end_run
    end if

    ! initialization

    ! modules warm up
    call common_info_init(namelist_file)
    call space_manager_init(namelist_file)
    call time_manager_init(namelist_file)
    call dynamic_quantities_init(namelist_file)
    call mesh_manager_init(namelist_file)
    call dist_recon_init(namelist_file)
    call diagnostics_init(namelist_file)
    call sample_data_init(num_dim, boundary)

    ! read in initial data
    call dynamic_quantities_input

    ! initialize shallow water related physical quantities
    allocate(height_gradient(num_dim,num_parcel))
    allocate(divergence(num_dim,num_parcel))

    ! height field
    h_index = dynamic_quantities_get_index(height_name)
    allocate(h_tendency(num_parcel))
    call dynamic_quantities_get_pointer(h_index, h)

    ! wind field
    allocate(u_index(num_dim))
    allocate(u_tendency(num_dim))
    do i = 1, num_dim
        u_index(i) = dynamic_quantities_get_index(wind_name(i))
        allocate(u_tendency(i)%value(num_parcel))
    end do

    ! time integration

    do while(time_manager_continue())
        call dist_recon_approximate(h_index)
        do i = 1, num_dim
            call dist_recon_approximate(u_index(i))
        end do

        ! debug
        call mesh_manager_parcel_to_grid

        do i = 1, num_dim
            call dist_recon_derivative(h_index, i, 1, height_gradient(i,:))
            call dist_recon_derivative(u_index(i), i, 1, divergence(i,:))
        end do
 
        ! calculate tendencies
        h_tendency = 0.0
        do i = 1, num_dim
            h_tendency = h_tendency+divergence(i,:)
        end do
        h_tendency = -h*h_tendency

        do i = 1, num_dim
            u_tendency(i)%value = -gravit*height_gradient(i,:)
        end do

        ! check computational stability
!        ...
!        call time_manager_adjust_dt

        ! update position of parcel
        call dynamic_quantities_parcel_move(u_index)

        ! update dynamic quantities
        call dynamic_quantities_update(h_index, h_tendency)
        do i = 1, num_dim
            call dynamic_quantities_update(u_index(i), u_tendency(i)%value)
        end do

!        call diagnostics_main(h_index, u_index)

        call time_manager_advance

        call dynamic_quantities_output

        write(6, "(' time ', f)") cur_time
    end do

    call run_manager_end_run

end program shallow_water_main
